X>’ /? A/ 







0 'xn 

P-/ 



(N&3A-CC- 1 7 9t 88 ) TEE KAETIAN CLIMATE AND 
ENERGY fcALANCL MODE 15 Vi 1 H CC2/H2G 
ATMOSPHERES Final fief-ort, 1 Mar. 1984 - 
Mar. 198b (New York Or<iv. , New York.) 37 

C 3 C L 


3 1 
P 

038 G 3/S 1 


N87-1 1673 


J n c s 
4 3931 


f 


NEW YORK UNIVERSITY 

FACULTY OF ARTS AND SCIENCE 

DEPARTMENT OF APPLIED SCIENCE 



FINAL REPORT 

THE MARTIAN CLIMATE & ENERGY BALANCE MODELS 
WITH C0 2 /H 2 0 ATMOSPHERES 

Martin I. Hoffert 
Principal Investigator 
NYU/86-154 

NASA GRANT NAGW-573 

March 1, 1984 - March 31, 1986 

New York University 
Washington Square 
New York, N.Y. 10003 


f 


NEW YORK UNIVERSITY 

FACULTY OF ARTS AND SCIENCE 

DEPARTMENT OF APPLIED SCIENCE 



1 


SERSONRL CYCLES OF C0 2 RND R 1 MOSPHERIC CIRCULRTION ON MRRS: 

RN ENERGY BRLRNCE RPPRORCH 

by 

Martin I. Hoffert, Vijay K. Narayanan and UJarren H. Ziegler 


RBSTRRCT 

The analysis begins with a seasonal energy balance model (EBM) for 
Mars. This is used to compute surface temperature versus x = 
sindatitude] and time over the seasonal cycle. The core model also 
computes the evolving boundaries of the C0 2 icecaps, net 
sublimational/condensation rates, and the resulting seasonal pressure 
wave. Model results are compared with surface temperature and 
pressure history data at Viking lander sites, indicating fairly good 
agreement when meridional heat transport is represented by a thermal 
diffusion coefficient D ~ 0.015 UU m~ 2 K" 1 . Condensational wind 
distributions are also computed. 

Rn analytic model of Martian wind circulation is then proposed, as an 
extension of the EBM, which incorporates vertical wind profiles 
containing an x-dependent function evaluated by substitution in the 
equation defining the diffusion coefficient. This leads to a 
parameterization of D(x) and of the meridional circulation which 
recovers the high surface winds predicted by dynamic Mars 
atmosphere models (~ 10 m s -1 ). Peak diffusion coefficients, D ~ 0.6 UJ 
m~ 2 K~\ are found over strong Hadley zones — some 40 times larger 
than those of high-latitude baroclinic eddies. Ulihen the wind 
parameterization is used to find streamline patterns over Martian 
seasons, the resulting picture shows overturning hemispheric Hadley 
cells crossing the equator during solstices, and attaining peak 
intensities during the south summer dust storm season, while 
condensational winds are most important near the polar caps. 
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INTRODUCTION 

Data from Viking Mars orbiters and landers continues to provide a 
picture of a sister planet marked by profound seasonal variations in 
surface temperature, volatile gas exchange and occasional dust 
storms (Tillman, 1985]. Some of the more dramatic events of Mars 
seasonal cycle are evoked in this descripitive passage (Haberle, 1986a): 

Fall has come to the northern hemisphere of Mars. In the mid-latitudes the 
mean temperature is dropping past - 70 degrees Celsius; at the north pole it has 
already reached - 123 degrees, cold enough to freeze carbon dioxide, the 
principal constituent of the thin Martian atmosphere, ft cap of dry ice is now 
forming on the pole. Before winter ends ice will reach the 50th parallel. 
Meanwhile in the southern hemisphere winter has ended, and carbon dioxide is 
evaporating at the south pole. Along the edge of the retreating polar cap a 
sharp temperature contrast between ice and soil warmed by the spring sun is 
giving rise to strong winds. During the short but hot southern summer, as Mars 
makes its closest approach to the sun, the winds will lift dust off the surface in 
great swirling storms. Ultimately dust may envelop the entire planet. 

Much of the prior Viking data analysis, as uuelt a current planning of 
future Mars science missions, is concerned with understanding 
meteorological seasonal cycles from observations by spacecraft 
orbiters and landers. Important related questions are whether and 
how such cycles have varied over Mars geological history, and the 
role they play in volatile gas exchange. 

Models at different complexity levels are instructive in such studies. 
For example, energy balance models (EBMs) with atmospheric transport 
neglected, or schematically depicted by thermal diffusion coefficents, 
have proven useful in understanding surface thermal physics (Leighton 
and Murray, 1966; Murray and Malin, 1973; Kieffer et al„ 1977; ToonM. 
al.. 1980; Hoffert et al. . 1981; James and North, 1982; James, 1985), 
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whereas more computer intensive three-dimensional general 
circulation models (GCMs) were historically used to model seasonal 
wind patterns. (Pollack et.. al. 1976, 1981; Haberle et al. . 1979). R 
simplification based on the nearly zonally-symmetric Mars circulation 
was exploited by Haberle et al. (1982), who used a two-dimensional 
dynamic circulation model, including dust storm effects; but extensive 
computations were still required. 

In this paper, simplified dynamics is taken a step further. Employing an 
integral method with assumed vertical atmospheric profiles, only the 
numerical solution of the insolation-forced, one-dimensional unsteady 
heat conduction equation is required. This allows circulation patterns 
for different seasons and orbital parameters to be rapidly assessed. 

The analysis begins with an energy balance approach to derive a core 
numerical model for surface temperature and carbon dioxide evolution 
over Mars seasonal cycle. Seasonally-varying temperature and 
pressure wave data from Viking landers are compared with the core 
model predictions to establish the realism of this part of the model. 
Next, new parameterizations of thermally driven and condensation 
winds, and baroclinic eddy transport effects, are developed. R test of 
our simple circulation model's surface velocity against existing GCM 
runs at north summer is presented. Finally, our wind model is employed 
to estimate the circulation streamline patterns over the four seasons 
on Mars. 


SURFRCE ENERGV BRLRNCE 


The core seasonal cycle computer model for T s (x,t), where x is the 
sine of the latitude and t is the time, is conceptually similar to that of 
James and North (1982): It is derived from a zonally-averagea energy 
balance in the spirit of the original Leighton and Murray (1966) model 
for seasonal volatile evolution on Mars, but includes the effect of 
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horizontal atmospheric transport on this balance. The core program 
yields as by-products the time-varying boundaries of north and south 
CC >2 icecaps and the rates of mass fluxing to and from the atmosphere 
from phase changes over the seasonal caps. Model parameters 
adopted here are summarized in Table 1. 

The basic energy balance involves four types of zonally-averaged 
vertical energy flux per unit surface area in a ground/atmosphere 
column: 


H sto^ x '^ + H subf x »^ = H rad( x >t) + ^hor^ x, ^» H) 

where H^o^t) is an energy sink by column heat storage, H su p(x,t) is 
an energy sink to sublimation of an underlying CO 2 frost, H rac j(x,t) is a 
net (solar minus longwave) radiant energy source and Hp or (x,t) is an 
energy source from horizontal atmospheric transport. 


Thermal Storage and Sublimation Heat Sinks 

Heat is stored during the evolving seasonal cycle by the thermal 
inertia of the atmosphere and land. The atmosphere mixes fast enough 
relative to seasonal timescales for the entire column heat capacity to 
come into play. Its heat capacity per unit area is therefore C a ^"p s c f? 
« 0.16 MJ m -2 K -1 , from T able 1 values. 

The effective heat capacity of the underlying ground depends on how 
deeply the seasonal wave penetrates, with pg, Cg, and A the ground 
density, specific heat per unit mass, and thermal conductivity, an 
annual heating wave of circular frequency 0 = 2rt/x p « 1.06 x 10~ 7 
rad-s _1 produces a thermal wave penetrating to depth d ^ 
[2A/(pgCgQ] 1/2 . The corresponding ground heat capacity per unit area 
is (Hoffert and Storch, 1979) 

C 9 * Pg c g d = 2A/(Od) = (2/ Q) 1/2 I, 
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where I = ipgCgAj 1 ^ 2 = (0/2j" /2 Cg is the ground thermal inertia. Notice 
that Cg scales with fT 1//2 , whereas I is a property of surface materials 
only. A typical ground thermal inertia on Mars is (Pollack et al .. 1981) I ^ 
300 J m -2 K -1 s -1//2 , corresponding a heat capacity per unit area of Cg ^ 
1.30 MJ m“ 2 K _1 for the seasonal wave. 

The total energy sink from column heat storage is therefore 

H st0 (x,t) ^ C9T s (x,t)/at, (2) 

where C = C a + Cg ^ 1.46 MJ m -2 K - ^ is the heat capacity per unit area of 
the atmosphere plus the thermally interactive ground column. 

A unique feature of Mars surface physics over solid CO 2 icecaps is 
sublimation to the gas phase of the main atmospheric constituent when 
T s > T c (p s ). Conversely, the atmosphere condenses to form C0 2 
icecaps when T s < T c (p s ). The pressure-dependence of the frost point 
temperature is approximately T c (p s ) 2 ; - ©|_/ln(p s /p*) ( where p* ^ 9.74 
x 10 11 Pa is a pressure scale, and ©|_ s; 3148 K is a temperature scale 
characterizing the latent heat of sublimation: L = R6 L ^ 5.95 x 10 5 J kg -1 
(James and North, 1982). Because of the weak logarithmic dependence, 
as surface pressure p s varies seasonally over its annual range 
(600-800 Pa, at the reference ellipsoid elevation), T c varies over a 
relatively much smaller range (149-151 K). For simplicity, a constant T c 
« 150 K is assumed over carbon dioxide icecaps. The heat sink to 
sublimation phase change over C0 2 icecaps is 


H sub^'^ - Lm(x,t). (3) 

where m(x,t) is the mass per unit time and per unit area injected to 
(removed from) the atmosphere by phase changes at the surface. This 
term is nonzero only over C0 2 ice when T s <. T c . 

Radiation Forcing 

Seasonal climate cycles are driven by an insolation function S(4>, L s ), 
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where <j> is the latitude and L s the heliocentric iongituae of the sun — 
the angular displacement along Mars' elliptic orbit measured from 
northern hemisphere spring equinox, L s □ 0. The eccenticity e, obliquity 
i and heliocentric longitude of perihelion L S p are parameters of the 
seasonal cycle, and are available for present epoch from Table 1. The 
diurnally- averaged incident solar radiation per unit area is given by 
(Berger, 1978): 

S(4>,L s ) ^ {S 0 [1 + ecos(L s -L S p)] 2 (Hsin(i)Sin£ + cos<i>cos<5sinH)}/[TT(1-e 2 ) 2 ], 

where d[l s ) is the solar declination angle given by sinS = sini sinL s and 
H(4>,L S ) is the hour angle for sunrise given by cosH = - tan4>tan<5; Here H = 
0 and H = rr correspond to polar night and polar day conditions 
respectively. 

Because Mars accelerates along its orbit near perihelion, and 
decelerates near apehelion, L s is nonlinearly related to elapsed time t 
from some reference point on the orbit. Neglecting terms of order e 2 , 
the time from northern spring equinox is related to L s approximately 
by 

t(L s ) ^ Tp/(2Tr){L s - 2e[sin(L s - L S p) + sinL S p]}. 


It is also convenient to work with a stretched latitude variable 
proportional to the area of the planetary surface from the equator, 

x(4>) = sin4>. 


These relations define S(x,t) for specified e, i and L S p. 

ft vertical column containing the atmosphere and thermally interactive 
ground is heated by a net radiative flux H raC j equal to the absorbed 
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solar flux less the longwave radiation flux lost to space. This term is 
modeled by 

Hradt*-^ = SIx.tKl - a) - ecr[T s [x,t]] 4 , (4) 

where cr ^ 5.67 x 10~ 8 W m~ 2 K“ 4 is the Stefan-Boltzmann constant and 



'^q > T s (x,t) > T c ' 


' £ g » T s^x-t) > T c ' 

Qr(X,t) = 

, e(x,t) = 


^i ; T s (x,t) <. T c , 


; T s (x,t) <T C , 


are planetary albedoes and emissivities, with different values over 
bare ground and ice. A CO 2 ice albedo of » 0.8 is assumed over both 
caps in the program. This is a rough average of north and south cap 
CO 2 albedos inferred from Viking data by Paige and Ingersoll (1985, Fig. 
53. A corresponding ice emissivity of Cj ^ 0.58 is needed to recover T c 
~ 150 K at the poles. 


Horizontal Heat Transport 

Horizontal (meridional) energy transport arises when winds blow fluid 
parcels containing enthalpy and gravitational potential energy across 
latitude circles. The meridional energy flux per unit area is pv(CpT + 
gz), where p, v and T are the density, meridional velocity and 
temperature; Cp is the constant-pressure specific heat, g the 
gravitational acceleration and z the altitude. The differential 
cross-sectional area perpendicular to the meridional energy flux can 
be written a(1 - x 2 ) 1//2 dzdA, where A is the longitude and x = sin4>. 
Accordingly, a general, three-dimensional, atmospheric circulation 
creates an energy flow per unit time across a latitude circle of: 

f2TT *• 


F(x,t) = a(1-x 2 ) 1/2 


OO 


pv(C D T + gz)dzdA. 


(5) 


Since a differential surface area element is dA s = 2TTa 2 dx, the net 
energy source per unit surface area into a column of surface area dA s 
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from horizontal flows is H^Q r = -aF/dR s = -(2TYa 2 j 'dF/dx; or 


H hor^ x ^^ 


inhere 


a 

9x 


aiy 

D(1-x 2 ) — 
ax 

/ 


- F(x,t) 

D(x,t) = 

2rra 2 (1 - x 2 )aT s /ax 


( 6 ) 


( 7 ) 


is the horizontal diffusivity coefficient for meridional energy flux. 

R constant D ~ 0.015 w m -2 K~1 independent of season or latitude was 
assumed in our initial numerical calculations. Rs discussed below, this 
value gave reasonably realistic surface temperatures and seasonal 
pressure wave distributions in comparison with Viking lander data. 
However, equations (5] and (7] are crucial links between any known 
circulation field and the thermal diffusivity. They will be employed 
subsequently to parametrize D[x] in terms of Mars wind circulation 
fields. 


Numerical Solution 

Rt constant D, substituting equations (2], (3), (4) and (6] into equation (1] 
gives an unsteady "heat conduction on a sphere" partial differential 
equation: 


9T< 


at 


i a 
x ax 


ai< 


(1-x 2 ) 


ax 


S(x,t)(1 - <*) - £cr[T s (x,tj]4 |_ni[x,t) 


( 8 ) 


i • > L-» a i / r~\ 

uui ioi e X = u/ u 


A /“S iIaP _ 

l.u X IU~ b 


i.D'C p is a timescale for horizontal 


diffusive heat transport. Numerical values are for D - 0.015 UJ m -2 K -1 . 
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Equation (8) is solved numerically for T s (x,t) by a Crank-Nicholson 
finite-difference scheme with 9T s /dx and 9 2 T s /3x 2 approximated by 
centered differences at interior points, and one-side differences at the 
poles, ujhere all space derivatives are evaluated at mid-timestep. The 
implicit Crank-Nicholson PDE solver employs computationally 
economical tri-diagonal matrices, is unconditionally stable, and has a 
second-order in space accuracy depending on (3 = At/KAx 2 ) (Ferziger, 
1981, p. 148). Normally, we want (3 <. 1, corresponding to timesteps At < 
(tAx 2 ); that is, the number of timesteps per Mars year needed for a 
given level of accuracy increases with the square of the horizontal 
resolution, we use a grid of 40 points in x = sin<t> (Ax = 0.05) as a 
compromise between small Ax desired to resolve icecap boundaries, 
and large Ax desired to minimize timesteps/year, and thus computing 
time, at given accuracy. To insure accuracy 1000 steps/year (At 
=1x10 -3 Tp) are used corresponding to (3 ~ 0.2. R D-coefficient 50 times 
larger would give (3 ~ 10. Some 10,000 steps/year at the same Ax are 
would be necessary to get (3 - 1 for an accurate solution. 

The program is initialized with a uniform annual mean temperature 
distribution and "spun up" for several Mars years to eliminate starting 
transients. It employs the following procedure to compute icecap 
boundaries and mass fluxes from sublimation/condensation over the 
caps: Normally Tim = 0 in equation (6) whenever the computed T s (x,t) > 
T c . when the computed surface temperature droped below the frost 
point, the program sets T s (x,t) = T c . The time-dependent cap 
boundaries x^(t) and xg(t) are evaluated as the gridpoints in the 
northern and southern hemisphere, respectively, where equation (7) 
predicts T s = T c . Since space and time derivatives of temperature are 
small over the constant-temperature caps, the sublimation mass 
source is simply: 

m(x,t) ^ {S(x,t)(1 - cq) - ejcrT c 4 )/L. (9) 

Balancing net radiant energy against phase change energy is an 
approximation dating back to the Leighton and Murray (1966) model. It 
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does seem justified however in light of detailed studies of Mars CO 2 
cap energy budgets by Paige and Ingersoll (1985), who found small 
contributions by transport and storage. Rates of mass injection to the 
atmosphere by south and north cap sublimation are computed by 
integrating the mass fluxes per unit area over the respective cap 
areas, 


x s (t) +1 

Mg(t) = 2rra 2 J m(x,t)dx, M N (t) = 2rra 2 J m(x,t)dx, (10) 

-1 X N (t) 

where negative values indicate mass removal by condensation. 


SERSONRL MODEL RESULTS RND VIKING DRTR 

Figure 1 shows the seasonal variation of pole-to-pole surface 
temperature computed from the energy balance model for Table 1 
parameter values. The results are expressed as constant surface 
temperature contours at AT s = 5 K intervals on the (x, L s ) plane. Rlso 
shown is latitude as a nonlinear vertical scale, <t> = 4>(x), and time in 
sots from N. Spring equinox as a nonlinear horizontal scale, t = t(L s ). 

The hemispheric asymmetry over the seasonal cycle is a consequence 
of maximum orbital distance from the sun (apehelion) and minimum 
orbital distance (perihelion) occuring presently at L s ^ 70° and 250°, 
respectively. These are sufficently close to the southern winter and 
summer solstices to markedly reenforce the seasonal thermal contrast 
in the southern hemisphere and diminish it in the northern hemisphere. 

The coldest spots on the planet are the CO 2 icecaps at 150 K; the 
warmest is the southern hemisphere hot spot at 240 K, extending as 
far south as 4> « 50 °S during southern spring and summer. Figure 1 
shows the computed seasonal boundaries of the north and south 
carbon dioxide icecaps x^t) and xg(t), assumed to form at the 
sublimation frost point temperature, T c ^ 150 K (The "wiggles" in the 
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receeciing cap boundaries are artifacts of the finite Ax resolution of 
the model). The CO 2 caps are asymmetric with the southern cap 
extending more equatorward than the north during their respective 
winter maximums. 

The model predictions of seasonal temperature histories at the 
appropriate latitudes can be directly compared with observations of 
surface air temperature made by VL-1 and VL-2 landers during the 
Viking mission. The latitudes of VL-1 (22 °N — now the Thomas Mutch 
Memorial Station) and VL-2 (48 °N) are indicated in Figure 1. Viking 
lander barometric pressure data at these sites exhibits the well-known 
seasonal CO 2 pressure waves associated with icecap cycles. The 
existence of condensation/sublimation pressure waves on Mars was 
predicted a decade before Viking by Leighton and Murray (1966), who 
assumed an energy balance model with no atmospheric transport. 


Lander Temperatures 

Figure 2 compares Mars surface temperatures measured at VL-1 and 
VL-2 over the seasonal cycle (Seiff, 1982) with predictions of the 
present energy balance model at the lander latitudes. 

The rapid temperature jumps computed during northern spring and 
summer near the VL-2 site are generated by abrupt sublimations of the 
receeding north CO 2 polar cap over one spatial step (Ax ^ 0.05), which 
propagate south by atmospheric diffusion (Cf. Figure 1). These appear 
to be artifacts of the finite-difference resolution, and converge to a 
smooth curve as Ax -*• 0. With this minor exception, the agreement 
between the model and the lander observations in both phase and 
amplitude for the seasonal thermal waves is remarkably good at both 
sites. 


Lander Pressures 

A more formidable challenge to the model is its ability to reproduce 



12 


the condensation/sublimation seasonal pressure wave observed by 
the Viking landers. Tillman (1985) has observed that interannual 
variations in the pressure wave at a given lander location are smalt. 
However, even after filtering out short-period oscillations associated 
with baroclinic "weather" disturbances, the VL-2 pressure wave of any 
given year is systematically higher than that at VL-1 (Seiff, 1982). This is 
mainly due to topographic differences at the two sites (VL-1 and VL-2 
are at elevations z\ - -1.5 km and Z 2 = -2.5 km, respectively, relative to 
the Mars reference ellipsoid). 


UUe therefore corrected the lander pressures pj(t) in Seiff (1982) to 
Mars’ reference ellipsoid elevation by using the hydrostatic equation: 
p s j(t) pj(t)exp( + Zj/hj(t)], where j = 1, 2. The effects of variations in 

scale height hj(t) = RT s j(t)/g from different surface temperature 
histories T s j(t) at the two sites were included in the correction. Since 
the Zj's are negative, surface pressures are lower than lander 
pressures in both cases, with a larger pressure correction for VL-2 (in a 
deeper valley). The result, shown in Figure 3, is that the VL-1 and VL-2 
pressure curves become virtually indistinguishable when these 
corrections are made. Rlso shown as an alternate vertical scale in 
Figure 3 is the atmospheric column mass loading at the "surface" 
reference ellipsoid elevation, m = p s /g. 


To model the seasonal pressure wave, we assumed mass added by 
phase changes over both caps from equations (10) is deposited in the 
atmosphere between the cap boundaries. Pressures over the caps are 
assumed to be maintained constant at the condensation frost point 
temperature T c . The rate of change of atmospheric mass between the 
caps is then 


M T (t) h 2TTa 2 (x N - x s ) 


a 

at 


»oo 

pdz = M^t) + 
JO 



( 11 ) 


independent of x. For the special case where M^ = - Mg the net flux 
into the system is zero; otherwise there is an increase or decrease of 
atmospheric mass. The surface pressure is assumed to adjust 
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"instantaneously" planetuuide to mass sources, and is found by 
numerical integration of 

dp s 9M T (t) 

— = ( 12 ) 

dt: (2TTa 2 ][><f N |(t) - Xg(t)] 


The results of this calculation are shown in Figure 3, where they can 
be directly compared with the lander-observed pressure wave 
corrected to the reference surface elevation. The pressure peaks 
result from more rapid degassing by the subliming pole than gas 
absorption by the condensing pole. The stronger peak in southern 
summer a consequence of hemispherically asymmetric heating 
augmented by atmospheric transport. Indeed, a feature supporting the 
adopted D - 0.015 UU m~ 2 K~\ at least at high latitudes, is its recovery 
of the observed higher amplitude south summer pressure peak. 


James and North (1982) did not obtain the observed relative amplitudes 
of the seasonal pressure peaks with their "clear sky" diffususion 
coefficient D ~ 0.004 w m -2 K _1 , although acceptable agreement was 
found when diffusivity was increased in their model to values as high 
as D - 0.1 UJ rrf 2 K _1 during 1977 global dust storms. Their inference is 
that D correlates with atmospheric opacity, increasing by orders of 
magnitude during large optical depth global dust storms. Such storms 
occurred in south summer of the first Viking year (Figure 2), precisely 
when high atmospheric heat transport is needed to rapidly sublime the 
south cap. But their correlation with optical depth may be fortuitous, 
since subsequent years when global dust storms did not recur exhibit 
essentially the same relative winter and summer peaks in the pressure 
wave (Tillman, 1985). By chosing a constant D ~ 0.015 UJ m -2 K" 1 — 
between James and North's "clear sky" and "optically thick" values — 
we get a fairly realistic pressure curve without invoking the dust 
storm connection. Even better agreement with the data, based on 
tuning the relevant paramers within observational constraints, seems 
possible, but was not pursued. 
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While the surface temperature curves at Viking lander latitudes are 
well represented, they may be too insensitive to transport effects to 
calibrate D. The seasonal pressure wave on the other hand seems 
mainly influenced by heat transport in the immediate vicinity of the cap 
boundaries. Clearly, it would be illuminating to if the thermal diffusion 
coefficient over the planet as a whole could be modeled in terms of 
dynamic atmospheric processes. This is the subject of the balance of 
this paper. 


MERIDIONAL CIRCULATION 

Unlike the earth's atmosphere for which extensive observational 
statistics on wind patterns are available (Oort, 1983], Mars atmospheric 
circulation is presently understood primarily from a few GCM 
simulations over restricted L s ranges. Existing lander anemometer 
measurements (Tillman, 1985], surface wind streak patterns (Kahn, 
1983], and dust storm observations (Haberle, 1986b] provide only 
tantalizing hints of the realism of these simulations. At this point 
therefore, simplified, partly analytical, models are useful tools for 
gaining insight into atmospheric circulation processes. 

Mean and Eddy Flows 

A reasonably good approximation is that the time for the atmosphere 
to radiatively equilibrate with the seasonally evolving surface 
temperature is sufficient short — about 10 sols according to Pollack 
OL (1981] — such that the circulation during the seasonal cycle of 
duration 670 sols is quasi-steady. The time, t, is then a parameter, so 
time-dependences of flow variables are omitted below for simplicity. 

It is also useful to distinguish between zonally-symmetric components 
of the circulation; and those which vary explicitly with the longitude, 
A. Any flow variable can be decomposed into such zonally-symmetric 
and longitudinally-varying components. Velocity and temperature, for 


15 


example, can be written: v(x,z) + vix,z,aj and T(x,z) + nx,z,A), where 
primes are departures from the zonal mean. Zonal averages around 
latitude circles of v' and T’ individually vanish identically; but the v'T 
product average (the meridional eddy flux) does not. On Mars, it is 
useful to further decompose v(x,z) into components from 
condensational winds, and thermally-driven cellular overturning in the 
meridional plane, 

v(x,z) = v c (x) + v T (x,z). (13) 

For simplicity, a vertically uniform distribution of condensational wind 
v c (x) is assumed. The thermally-driven component vy(x,z), reverses 
direction at higher altitudes forming closed Hadley-like cells, and is 
therefore altitude-dependent. It has the property Jpv ydz = 0. Both 
components transport energy across latitude circles, but only the 
condensational wind transports mass. 

Let meridional diffusion effects from zonally-symmetric processes be 
represented by a mean-flow diffusion coefficient D m , those from 
eddies or waves in the east-west direction by an eddy diffusivity D e . 
(Logically, D m too can be considered an "eddy diffusivity" associated 
with the mean flow overturning eddies of the Hadley circulation in the 
x,z plane, but no confusion should result if the terms are understood.) 
It follows from equations (5) and (7) that the diffusion coefficient is 
expressible 

D(x) = D m (x) + Dg(x), (1*4) 

where 


oo 

- Jpv(CpT + gz)dz 
o 

i 

a(1 - x2j)/2 aT s /0x 


D m (x) = 


(15a) 
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2rr OO 

- I J pCp(v'r)dzd^ 

D e (x) ) = — 0_Q . (15b) 

2ira(i - x 2 ) 1/2 aT s /ax 


The eddy component from unstable baroclinic waves can be estimated 
from Stone's (1972) baroclinic instability parameterization. Its value for 
the planet as a whole is approximately (Hoffert et al. . 1980) 

0 . 28 R 2 p s CpT s ,/ 2 (r ad - n '/2 (I s - T c ) 

D e ~ 0.018 W m~ 2 K" 1 , 

-n-g3/2 a 3f2 

where the numerical estimate corresponds to Table 1 parameter 
values and a Coriolis parameter f « 1.4 x IQ" 4 s~1 jt j S close to the 
diffusion coefficient value used in the core model (D - 0.015 W m" 2 K" 1 ), 
and may characteristizes transport near the cap boundaries where 
baroclinic activity dominates (Hablerle, 1986b). This may be why the 
seasonal pressure wave was modeled fairly well. 


In contrast to the terrestrial circulation, Martian winds are apparently 
dominated at low latitudes by the zonally-symmetric circulation 
(Haberle et al.. 1982). UUe focus next on these components of the 
circulation, for which no dynamically-linked diffusion parametrization 
has been developed as yet. 


C.Qnagn£3tiQ na l w i nd 

The mass flowing across a latitude circle carried by the condensational 
wind is: 


M c (x) = 2rra(1-x 2 ) 1//2 


>oo 

pv c dz. 


(16) 


JO 

This quantity appears in the unsteady continuity equation integrated 
from the surface to the top of the atmosphere: 



17 


a 


POO 


9tJ 


pdz 

0 


i a 
a ax 


✓ 

(1-x2)1/2 

"OO 

pv c dz 

k 0 

0 


= 0 . 


(17) 


Substituting equations (11) and (16) into (17) yields an expression for the 
x-wise derivative of the condensational flow, 6M c /6x = M y/(xg - x^). 
UUe therefore use the linear distribution, 


M c (x) = (M N + Mg) 


'x - x N ' 

*s -xn, 



(18) 


which has the proper derivative, and satisfies the cap boundary 
conditions M c (x N ) = - M N and M c (x s ) = Mg. 


The condensational velocity distribution is obtained by integrating 
equation (16) with the density profile p(z) ^ p s e _Z//h , where h = p s /p s g 
is the scale height: 


9M c (x) 

v c (x)= 

2rrap s (1-x 2 ) 1 / 2 


(19) 


Condensational winds versus latitude at the equinoxes and solstices 
computed from equations (18) and (19) using the sources M N and Mg 
from the core program are illustrated in Figure 4. The peak 
condensational wind magnitude, v c _ 0.4 m s~\ occurs near the south 
cap boundary in northern winter (L s = 270°) as this cap sublimes under 
the influence of high southern hemisphere insolation and heat transfer 
to the cap boundary. This wind is comparable in magnitude and location 
to peak values computed by James (1985, Fig. 3). Rt this L s , the 
condensational flow is everywhere northward, depositing mass in the 
atmosphere along the way, and contributing to the northern winter 
peak of Fig. 3, until what remains condenses at the north polar cap. Rn 
opposite but weaker flow exists in northern summer (L s = 90°), when 
southward flow generated by north cap sublimation is absorbed by the 
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south cap. if sublimation and condensation rates mere precisely equal 
at opposing caps, condensational winds uuould still blow, but seasonal 
pressure waves from this effect would be absent. These winds are 
weaker at the equinoxes, as expected, but at L s = 180°, v c reverses 
direction mid-planet to accomadate the model prediction of 
condensation at both caps during this time. 

The fact that average condensational winds are of order v c 0.2 m s _1 
raises an interesting question. The transit time for a gas parcel 
sublimed at one pole to cross the planet to is rra/vc ~ 5.3x10 7 s - 600 
sols, so the time to physically transport condensational mass between 
VL-1 and VL-2 lander locations is [[48-22)/180] x 600 sols ~ 87 sols. Why 
then is there no phase lag of this order in the pressure waves 
observed the two stations? Recall that both lander pressure curves, 
corrected to the reference ellipsoid and with baroclinic weather 
oscillations removed, are nearly coincident. One possibility is that the 
atmosphere rapidly adjusts on the diurnal timescale of internal 
gravity waves to create horizontally uniform surface pressures. The 
observational absense of a discernable phase lag seems in any event 
a sufficient rationale for the instantaneous pressure adjustment 
assumed in equation (12). 


Thermally-driven winds 

The zonally-symmetric thermally-driven circulation is composed of a 
mean merdional circulation, in which the merdional v(x,z) and vertical 
w(x,z) velocities are coupled by continuity, and a zonal mean wind 
u(x,z) and temperature T(x,z), which are coupled through the thermal 
wind relation. These components are in turn coupled through the 
Coriolis torque. Advection of relative angular momentum in inviscid 
models is also necessary to describe thermal forcing of the mean 
meridional circulation. 

For nearly inviscid models, the Hadley cell is limited in its latitudinal 
domain, which in turn limits the intensity of the zonal flow. A major 
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difference between terrestrial and Martian zonal mean flows is the 
the relative absense of eddy momentum fluxes on Mars, which makes 
the circulation nearly zonally symmetric over ± 30° latitudes (Leovy 
et at. . 1973]. For example, a two-dimensional zonally symmetric 
circulation model, such as that employed by Haberle et al. . 1982] to 
study Mars dust storms, would probably not describe the earth's 
atmospheric circulation very well, fis a consequence of the earth's 
land-sea contrasts eddy momentum fluxes are important at all 
latitudes, and need to be specified from observations in 
zonally-averaged models to recover a realistic circulation 
(MacCracken and Ghan, 1986]. 

To model D m an integral method was used in which a meridional wind 
shear of the form 9v/9z ^ - v(x]/h is assumed. This is formally 
analogous to the zonal wind shear of the thermal wind equation, 
9u/9z ss u(x)/h, where u(x] = (gh/afT s ]9T s /94> (Washington and 
Parkinson, 1986]. Rt this point v(x] is an unknown function to be 
determined in terms of D m (x]. However, it is well-known that latitudinal 
heating gradients drive the zonally-symmetric (u,v] circulation in 
general (Holton, 1979; Haberle, 1986b]. In particular, it follows from 
equation (15] of Leovy et al. (1973] that for steady, inviscid flows in 
equatorial regions, and neglecting eddy momentum fluxes, v(x] ~ - 
u(x). UJe therefore assume as a crude approximation that 

v(x] ~ - (gh/afT s )(1 - x 2 ] 1/2 9T s /9x (20] 

in latitude bands of the Hadley circulation. 

The linear profile v-p(x,z] = - v(x][1 - (z/h)] also satisfies zero mass flux 
across latitude circles, Jpv-j-dz =0 for p ^ p s e _z / h , and produces a 
thermally driven surface wind flowing toward local hot spots (v s = 
vj(x,0) oc 0T s /9x). This wind reverses direction above the scale height 
corresponding to overturning Hadley cells, but needs to be carefully 
interpreted at high altitudes where unrealistically large negative 
velocities are implied as z -*• oo. The linear profile seems a reasonable 
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first approximation in the lower troposphere (z < zj), where Mars' 
"tropopause" is at Zj « 40 km. The temperature at zj is Tj ^ 160 K. 

The scale height is expressible in either of the alternate forms 
RT s (xJ p s 

h(x) « = , (21) 

9 Ps^9 

where typical planetary mean value is "h » 11 km. The troposphere 
lapse rate assuming linear temperature profiles is 
9T T s (x) - T t 

r(x) = « . (22) 

9Z Zj 

The planetary mean lapse rate in Mars' stable atmosphere is T « 1.45 
x 10-3 K rrfl; whereas the adiabatic lapse rate is r ac | = g/Cp « 4.55 x 
10~3 K m _1 > T(x). It is useful to define a dimensionless relative stability 
parameter, 


r(x) = (R/g)[r acj - r(x)] 

(23) 

To summarize, the vertical structure 
approximated by the profiles, 

of the atmosphere is 

p(x,z) « p s (x)exp{-z/h(x)}, 

(24) 

T(x,z) « T 0 (x) - r(x)z, 

(25) 

v(x, z) « v c (x)-v(x){1-[z/h(x))}. 

(26) 


Substituting equations (21H26) into equation (15a) and performing the 
required integrations yields: 


v(x) = 


fga(1-x 2 ) 1|/2 9T s /9x'| 


v c (x)(1 + rtx)} 
Yix) 


l p s CpT s (x)y(x) 


* D mM 


(27) 
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UJe have departed in equation (27) from the Leovy et al. (1973) integral 
method in two ways: (i) v(x) is evaluated in terms of D m (x) as opposed 
to atmospherically-absorbed radiation and (ii) the condensational wind 
effect is included. 

Neglecting condensational winds, and eliminating v(x) between 
equations (20) and (27) with the help of equation (21), gives the 
following estimate of the zonally-symmetric thermal diffusion 
coefficient of Hadley cells: 

~P s c p9^ _ 

D m (h/a) 2 - 0.6 W m -2 K -1 , 

f 

using Table 1 values. This is more than 30 times larger than the 
baroclinic eddy value, D e , estimated previously. Of course, we would 
not expect this to apply beyond the x's bounding central Hadley cells, 
so a rather strong drop-off to values characterizing the cap region is 
implied. 

Both two and three-dimensional primitive equation general circulation 
models indicate that during solstice conditions cross-equatorial Hadley 
cells exist surface winds of order v s - 10 m s -1 over latitudes within + 
30° of the equator (Haberle et al.. 1982). Typical surface wind 
distributions, v s (x), during south summer computed with 
two-dimensional (2D, Haberle et al.. 1982) and three-dimensional (3D, 
Pollack et al.. 1976) dynamic models are shown in Figure 5. These 
dynamically explicit models are identical mathematically and 
numerically, except that ^-derivatives are set equal to zero in the 2D 
model, automatically suppressing any longitudinal waves or eddies. 
The closeness of the 2D and 3D results supports the assumed 
dominance of the zonally symmetric winds. These zonal mean 
meridional surface winds on Mars are substantially more intense than 
their terrestrial counterparts — the latter being of order v s ~ 1 m s _1 
(Oort, 1983, p. 157). 
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UUe sought to determine whether a D m (x) distribution could be found 
with peak values of D ma>< ~ 0.6 UJ rrr 2 k _ 1 which approximates the 
surface wind distributions of the dynamic models. To test this we used 
T s , 9T s /9x, y and v c versus x from the core model at L s = 90° — close, 
but not identical, to the GCM L s « 106° conditions — and evaluated our 
model's corresponding surface velocity 


v s (x] = 


rga(l-x 2 ) 1/2 9T s /ax'i 


l Ps^pT (x) J 


D m [x) 


v c (x]{1 + 2r(x)} 


rix) 


for various assumed functions D m (x). 


fin immediate finding was that surface wind distributions consistent 
with GCM results at all latitudes could not be obtained with a 
latitude-independent diffusion coefficient. R constant D m , matching v s 
~ 10 m s _1 over Hadley cell regions gave too large surface winds near 
the cap boundaries (where 9T s /9x is large], while a small D m could 
recover cap region surface winds, but not those of Hadley zones. 
Clearly, a latitudinally-varying diffusion coefficent is required. 


By experimenting with different functions, a reasonably successful 
diffusion coefficient distribution was found. Rs shown in Figure 5, the 
model produced a v s (x) distribution close to those of the dynamic 
models — almost as close as the 2D and 3D surface wind distributions 
are to each other — when the logarithm of D m varied parabolically 
with x between the cap boundaries. This was a very welcome finding, 
as it was not even certain a priori that a diffusivity approximation 
would appropriately describe the zonally-symmetric circulation. Rs it 
happened, we were able to take the results even a further step to 
compute the streamlines. 


In Figure 5 [TOP] the D m (x) used can be compared with James and 
North’s "clear sky" and "dust storm" diffusion coefficients discussed 
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previously, and typical values of the Daroclinic eddy component D e . 
This D m (x) can be written 

lnD m (x) = lnD min [1 - 4Kf(x)], (28] 

where K = ((lnD ma>< /lnD mjn ] - 1] and f(x] = [(x - xg)(x - Xn)/(x n - x s ] 2 ]. 
Diffusivity is a minimum (D m j n = 0.015 UU m -2 K _1 ) at the cap boundaries 
Xfs| and xg, and peaks (D ma>< = 0.60 UU rrf 2 K -1 ) midway between the cap 
boundaries. Near cap boundaries, D m (x] is sufficiently small that 
surface velocity is dominated by condensational winds, v s ~ 
v c U + 2y]/Y- For y ~ 0.15, and v c ~ - 0.2 m s" 1 , this gives a southward 
surface velocity near the cap boundaries of v s - - 2 m s~\ close to 
values shown of Figure 5. Notice that the surface winds induced by the 
condensational winds exceed them by almost an order of magnitude, 
an effect produced in the dynamic models as well. 

Despite many approximations, the model reproduces GCM surface wind 
patterns near north summer solstice suprisingly well, including three 
overturning cells and condensation wind effects. Even details like the 
local surface velocity dip near the equator produced by the solar 
heating function are recovered. 

Rn interesting result was that the diffusion coefficient is much larger 
across Mars Hadley cells, than at the cap boundaries. R comparable, 
but latitude-independent D ~ 0.6 UU m~ 2 K~\ will match observed zonal 
mean surface temperature distributions on earth (Hoffert et al„ 1983) — 
also surprising in light of the earth's 100 times greater surface 
pressure. In the terrestrial case, this must result from very different 
flow mechanisms including atmospheric eddy transport and ocean 
currents. 


Streamlines 

In our flow model, equations (19) and (26H28) together with 
information from the core EBM, define the zonally-symmetric 
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meridional velocity field v(x,z) at any time during the seasonal cycle. 
These patterns are more clearly visualized by looking at the 
streamline contours, than at the velocity field as such. 


The partly analytic nature of the model makes streamline calculations 
fairly straightforward: A streamfunction ip(x,z) can be defined by its 
partial derivatives dfy/dx = -2rra 2 pw and dty/dz = 2rra(1 - x 2 ) 1//2 pv 
such that continuity is identically satisfied: 


a(pw) d [(1-x 2 ) 1/2 pv] 

+ = 1/(2rra 2 ) 

dz a 9x 


d 2 # 

azax 


a 2 # ' 
axaz 

✓ 


= o. 


The streamlines — contours of constant mass flow per unit time across 
which there is no flow — are therefore available from 


#(x,z) = 2rra[1 - x 2 ) 1//2 


rz 


p(x,z')v(x,z')dz'. 


0 


Substituting equations (24) and (26) and integrating gives: 


#(x,z) = 2rra(1 - x 2 ) 1/2 (p s /g){- v(xMz/h]-e" z/h + v c (x)*[1 - e“ z/h ]}. (29) 


Notice that the streamfunction at the top of the atmosphere equals 
the condensation mass flux ijj(x, «>) = M c (x), and that the surface and 
poles are zero mass flux boundaries. Tp(x,0) = ip(±1 ( z) = 0. 


Figure 6 shows the meridional plane streamlines computed at L s = 0°, 
90°, 180° and 270°. The contours of constant ijf have units of 
Megatonnes per second (1 Mt s" 1 = 10 9 kg s -1 ), and in the present sign 
convention are < 0 for clockwise circulations in the x,z plane. At north 
spring equinox we find north and south Hadley cells of roughly 
comparable intensity in each hemisphere with core flows centered at 
the scale height. A weak hemispheric asymmetry is produced even at 
equinoctal conditions because the ground thermal inertia produces 
asymmetric surface temperature distributions which "remember" their 
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history. 

As the northern hemisphere warms during its summer solstice, the 
south main cell migrates northward, penetrating to some 30° N 
latitude, and intensifies to flow rates of - 18 Mts' 1 ; fit the same time 
the northern cell weakens, and a weak indirect cell is produced near 
the south cap boundary by condensational winds. The Hadley cells 
become more nearly balanced in the two hemispheres at northern fall 
equinox, although again there is some hemispheric asymmetry. Finally, 
during northern winter solstice (southern summer), the northern cell's 
boundary moves southward, and attains a core circulation of - 28 Mt 
s -1 — the highest intensity of the seasonal cycle. Indeed, southern 
summer is the time of year on Mars when intense surface winds can 
attain the saltation velocities needed to raise particulate matter into 
the atmosphere. This is often a precursor of global-scale dust storms. 

CONCLUDING REMARKS 

UJe have seen that the core energy balance model developed in the 
early part of this paper gave a reasonably good account of the 
seasonal surface temperature and pressure waves on Mars as 
recorded by the Viking landers. Moreover, our extension to wind field 
analysis of the energy balance approach gives qualitatively plausible 
descriptions of atmospheric circulation over all seasons, and 
inherently recovers GCM-predicted surface winds at the north summer 
calibration point. The wind field depends partly on the pole-to-pole 
surface temperature distribution, and partly on the parameterized 
diffusion coefficient distribution. 

More work is needed on modeling the diffusion coefficient distribution 
in terms of limitations on the boundaries of the Hadley circulation from 
angular momentum and stability constraints. Moreover, a feedback on 
the core model is suggested by the wind field analysis: In particular, 
internally-consistent calculations of surface temperatures at 
equatorial latitudes should employ the D m (x) of the wind analysis in 
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the heat conduction equation, as opposed the the constant D ~ 0.015 UU 
m~ 2 K -1 assumed here. In turn, surface temperatures computed with 
this D m (x] feed back on the wind field. Prior to our wind analysis, we 
accepted the conventional wisdom of a largely decoupled lunar-like 
planet in which atmospheric transport feedbacks on surface 
temperature distributions are small and confined to polar cap 
boundaries. This assumption needs to be reexamined. 

In light of its computational economy compared to GCMs, applications 
of the present approach to related phenomena of the Mars seasonal 
cycle are promising areas of future research. For example, in the Mars 
seasonal hydrological cycle as modeled by James (1985], it is assumed 
that meridional transport of water vapor is accomplished by 
condensational winds plus some unspecified diffusionlike process. Our 
wind model suggests Hadley Cell transport, in which the smalt scale 
height of H 2 O profiles would selectively favor transport by 
near-surface winds, might provide an explicit explanation of this 
"diffusive" transport. Viking orbiter data on vertically-integrated water 
vapor are available for testing such models. 

Rnother potential application is modeling of Martian dust storms. The 
zonatly-symmetric circulations simulated by Haberle et al. (1982] at 
south summer conditions demonstrate that positive radiative 
feedbacks from dust raised by surface winds strengthens and expands 
the Hadley circulation. A puzzling question is why global dust storms 
are not triggered every Mars year. The answer may be related to the 
interplay of the zonalLy-symmetric and baroclinic eddy components of 
the circulation (Haberle, 1986b]. Although as presently configured our 
model does not incorporate the radiative feedbacks needed to mimic 
such global dust storms, such extensions can be implimented in 
principle. A simplified representation of zonally symmetric and 
baroclinic eddy transports, including the influence of vertical thermal 
structure, might well provide useful insights into the complex 
dynamics of dust storm events. 
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Finally, the present approach seems well suited to studies where runs 
over a large parameter space are needed. These include sensitivity 
studies, analyses of Mars paleoclimate when orbital parameters were 
different, and analysis of future scientific missions to Mars — where 
instrument design may require estimates of meteorological conditions 
over entire annual cycles. 
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T able 1. orbitrl rnd phvsicrl chrrrcteristics of mrrs (current epoch) 


Property 

Value 

Units 

Rnnual period, x^ = 2tt/0 

5.94 x 10 7 

(s) 

Diurnal period, T d , 1 sol 

8.86 x 10 4 

[S] 

Eccentricity, e 

0.0934 

[-] 

Obliquity (spin axis tilt), i 

25.1 

[deg] 

heliocentric longtude of perihelion, L S p 

250 

[deg] 

Solar constant, S 0 

586 

[uu m -2 ] 

Planetary radius, a 

3.39 x 10 6 

[m] 

Gravitational acceleration, g 

3.73 

[m s -2 ] 

Bare ground heat capacity per unit area, Cg 

1.30 x 10 6 

[Jm -2 K" 1 ] 

Bare ground thermal inertia, I 

300 

[Jm -2 K -1 s-1/2] 

Bare ground albedo, or g 

0.30 

[-] 

Bare ground emissivity, e g 

1.00 

[-] 

Ice albedo, o q 

0.80 

[-) 

Ice emissivity, £| 

0.58 

[-] 

Ice sublimation temperature, T c 

150 

[K] 

Latent heat of CO 2 ice sublimation, L 

5.95 x 10 5 

(Jkg- 1 ) 

Gas constant, R 

189 

(J kg -1 K" 1 ) 

Rtmosphere specific heat, c g 

830 

[Jkg -1 K -1 ] 

Mean surface pressure, "p s 

700 

[Pa] 

Mean mass per surface area, m =-~p s /g 

188 

[kg m -2 ] 

Mean surface temperature,T s 

218 

[K] 

Mean surface density, ~p s ="p s /RT s 

0.0170 

[kg m -3 ] 

Mean scale height,~h = rT s /g 

1.1 x 10 4 

[m] 

Tropopause height, z T 

4.0 x 10 4 

[m] 

Tropopause temperature, Ty 

160 

[K] 

Rdiabatlc lapse rate, r ad = g/Cp 

4.55 x 10“ 3 

[K m -1 ] 

Mean troposphenc lapse rate, T = (T s - T jVzj 

1.45 x 10" 3 

[K m -1 ] 

Mean stability parameter, y = (R/g)(r ad - T ) 

0.15 

[-] 
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FIGURE CRPTIONS 

Figure 1. Mars seasonal surface isotherms computed from the present 
model. 

Figure 2. Seasonal variation of near-surface atmospheric temperatures 
at Viking lander sites [Seif f , 1982) compared with present energy 
balance model predictions. 

Figure 3. Seasonal variation of atmospheric pressure (reduced to Mars 
reference ellipsoid) at Viking lander sites compared with present 
energy model predictions. 

Figure 4. Latitudinal distribution of condensational winds at north 
spring equinox (L s = 0°), north summer solstice (L s = 90°), north fall 
solstice (L s = 180°) and north winter solstice (L s = 270°). UUind directions 
are positive toward the North Pole. 

Figure 5. (TOP). Diffusion coefficient distribution D(x) used to match 
surface wind distribution of dynamic circulation models; (BOTTOM) 
Surface velocity distributions near north summer solstice. Present 
model results are for D m = D(x) at L s = 90°. 

Figure 6. Streamline patterns on Mars over the four seasons from the 
present model. 
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